On Holmboe's instability for smooth shear and density profiles 
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The linear stability of a stratified shear flow for smooth density profiles is studied. This work 
focuses on the nature of the stability boundaries of flows in which both Kelvin-Helmholtz and 
Holmboe instabilities are present. For a fixed Richardson number the unstable modes are confined 
to finite bands between a smallest and a largest marginally unstable wavenumber. The results in this 
paper indicate that the stability boundary for small wavenumbers is comprised of neutral modes with 
phase velocity equal to the maximum/minimum wind velocity whereas the other stability boundary, 
for large wavenumbers, is comprised of singular neutral modes with phase velocity in the range of 
the velocity shear. We show how these stability boundaries can be evaluated without solving for the 
growth rate over the entire parameter space as was previously done. The results indicate further 
that there is a new instability domain that has not been previously noted in the literature. The 
unstable modes, in this new instability domain, appear for larger values of the Richardson number 
and are related to the higher harmonics of the internal gravity wave spectrum. 

PACS numbers: 



I. INTRODUCTION 



Stratified shear flow instabilities occur in a variety of 
physical contexts such as astrophysics hi, the Earth's 
atmosphere and oceanography 0, HJ 0> HI • The linear 
instability problem of stratified shear flows has been ad- 
dressed in a large body of literature. Since the original 
work of Helmholtz and Lord Kelvin many models 
of flows and density stratifications have been investigated 
both analytically and numerically. Parameterized by the 
global Richardson number, these investigations have re- 
sulted in a large variety of stability/instability domains in 
the Richardson number - wavenumber space lESElEI- 
However, a full understanding of this large variety of sta- 
bility domains does not exist. 

A step towards understanding a particular aspect of 
shear flow instability was achieved by Holmboe [l^] • Us- 
ing a piece-wise-linear form of the velocity and density 
profile Holmboe managed to distinguish between two 
classes of unstable modes that are present in stratified 
shear flows. In the first class the unstable modes have 
zero phase velocity in a reference frame of zero-mean ve- 
locity and exist for a finite range of strength of the strati- 
fication (Richardson number). This instability is referred 
to in the literature as Kelvin-Helmholtz instability since 
the behavior of the unstable modes resembles the one pre- 
dicted by Kelvin-Helmholtz. The unstable modes of the 
second class have non-zero phase velocity and typically 
smaller growth rate than the Kelvin-Helmholtz modes. 
However they are present for arbitrarily large values of 
strength of the stratification, making them better candi- 
dates to explain certain physical phenomena. This sec- 
ond kind of instability is referred to as Holmboe instabil- 
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ity. 

Several authors have expanded Holmboe's work 0, 
ItH h~5| by considering different stratification and veloc- 
ity profiles that do not include the simplifying symme- 
tries Holmboe used in his model. Discussions on the 
mechanisms involved in the instability can be found in 
E]|^J. Hazel 01 an d more recently Smyth and Peltier 
[l8j have shown that Holmboe's results hold for smooth 
density and velocity profiles as long as the length scale 
of the density variation is sufficiently smaller than the 
length scale of the velocity variation. Further, effects 
of viscosity and diffusivity 0, |2(J , non-linear evolution 
pH I22I l2l| and mixing properties [24[ of the Holmboe 
instability have also been investigated. Experimentally 
Holmboe's instability has been investigated by various 
groups. Browand and Winant |25| first performed shear 
flow experiments in stratified environments under the 
conditions that Holmboe's instabilities are present. Their 
investigation has been extended further by more recent 
experiments [H mEHHH HI ■ The unstable Holm- 
boe modes have been observed and Holmboe's predic- 
tions verified. 

In this work we examine the stability of smoothly strat- 
ified shear flows as was done in [13 and 0, but here we 
focus on finding the kind of marginally unstable modes 
that comprise the stability boundaries. We show how 
these modes can be determined without solving the full 
eigenvalue problem for the complex eigenvalue c. Fur- 
thermore, our results indicate that new instability re- 
gions exist which have not been discovered before in the 
literature. 

This paper is structured as follows. In section [H] wc 
formulate the linear stability problem and discuss the 
possible marginally unstable modes that can comprise 
the stability boundaries. In section II I II we describe the 
numerical methods used. Section ITVI presents our results. 
Specifically, we show the location of the marginally unsta- 
ble modes and numerically verify that they indeed con- 
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stitute the stability boundaries. Conclusions are drawn 
in section [3 where we also give a physical description of 
our results. 



II. FORMULATION 

We begin with the Taylor- Goldstein equation, which 
describes linear normal modes of a parallel shear flow in 
a stratified, inviscid, non-diffusive, Boussinesq fluid: 



U" 



U-c (U-c) 2 



= 0, 



(1) 



where 4>(y) is the complex amplitude of the stream func- 
tion for a normal mode with real wavenumber k and 
complex phase velocity c. lm{c} > implies instabil- 
ity with growth rate given by 7 = fclmjc}. U{y) is 
the unperturbed velocity in the x direction and J(y) = 
—gp dp/dy is the squared Brunt- Vaisala frequency. 
Prime indicates differentiation with respect to y. Wc 
note that when c is real and in the range of U there is a 
height y c , at which U(y c ) — c. At the height y c , called the 
critical height, equation ||TJ has a regular singular point. 
Many of the features of the unstable modes are related to 
the presence of this singular point. Equation together 
with the boundary conditions 4> — > for y — > ±00, forms 
an eigenvalue problem for the complex eigenvalue c. 

The Taylor- Goldstein equation (JIJ, has been studied 
for many different density and velocity profiles and is 
known to have four different classes of modes as so- 
lutions [3l]]: (i) For some conditions unstable modes 
exist with the real part of the phase velocity within 
the range of U. The phase speed of these modes sat- 
isfies Howard's semi-circle theorem |c — l/2(sup{£/} + 
inf{[/})| < l/2|sup{[/}-inf{[/}|. Furthermore, if these 
modes exist the Miles-Howard theorem guarantees 
that somewhere in the flow the local Richardson number 
defined by: 



My) = 



[u'(y)Y- 



(2) 



must be smaller than 1/4. As discussed in the intro- 
duction, for some velocity and density profiles, the class 
of the unstable modes can be further divided into two 
subclasses of unstable modes, those whose phase veloc- 
ity is zero with respect the mean flow (Kelvin-Helmholtz 
modes) and those whose phase velocity is non-zero, the 
Holmboe modes, (ii) The complex conjugates of the un- 
stable modes constitute the second class of modes. That 
is, for each (c u , <f> u ) that is a solution of (JTJ that de- 
scribes an unstable mode there is a solution (q, <pd) that 
describes a damped mode with c<j = c* and 4>d — 4>* a - 
(iii) The third class of modes consists of non-singular 
traveling modes with phase velocity outside the range 
of U . These are internal gravity wave modes modified 
by the shear, (iv) Finally, due to the singularity at the 
critical height y c a continuum of singular neutral modes 



exist with a singular behavior of the first derivative of the 
steam function at the critical height. The phase veloc- 
ity of these modes lies within the range of U. For these 
modes an asymptotic analysis close to the critical height 
shows that the stream function can be written as the sum 
of two linearly independent solutions: 

<j>~A ± {\y-y D \*+...) + B ± (\y-y c \ s > + ...) (3) 



w here s a = 1/2 - ^1/4 - Ri(y c ), s b = 1/2 + 

— Ri(y c ) and the index ± corresponds to the so- 
lution above or below the critical height. The solutions 
in equation J3J are known as the Frobenius solutions. If 
these modes are marginally unstable the Frobenius solu- 
tions need to be considered as the limit Im{c} — > + and 
yield the connection formulas: 



A- = 



„ts a Tr A + 



A+ and B~ = 



(4) 



For each wavenumber and fixed density and velocity pro- 
files more than one of the above mentioned modes can ex- 
ist. Consequently, it is interesting to try and find bound- 
aries, if they exist, that separate these four classes of 
modes. This is what we try to achieve in what follows 
for a specific example of shear and stratification. 

The example we investigate follows Hazel's [l^ work 
where the velocity profile is given by: 



U(y) = tanh(y). 



(5) 



Note that the problem has been is nondimcnsionalized 
using the half-thickness and the half velocity change of 
the shear layer. For the stratification as, in we pick 
the squared Brunt- Vaisala frequency's functional form 
(in nondimensional units) to be given by: 



J(y) = J sech 2 (Ry) 



(6) 



where gives the non-dimensional length scale of the 
density stratification. The symmetry J(y) = J(—y) 
and U{y) — —U{—y) implies that if (c, (f){y)) is the 
eigenvalue and eigenfunction of an unstable mode, then 
(— c*, (j>(— y)) is also a solution of equation |QJ that de- 
scribes an unstable mode traveling in the opposite direc- 
tion. This allows us with no loss of generality to focus 
only on the modes that have Re{c] > 0. 

The following observations illuminate the origin and 
possible location of the stability boundaries. First, as 
noted by Hazel ^7], for the density and velocity profile 
given in equations (|6I5() respectively, the functional form 
of the local Richardson number crucially depends on the 
value of the parameter R (see figure [TJ. For R < \/2, the 
local Richardson number has a unique minimum at y = 
given by Ri(0) = Jo- This implies that we can only have 
instability if Jo < 1/4 |32| and any stability boundary 
must lie within this region or on its boundary. For \/2 < 
R < 2 two minima exist symmetrically around zero and 
therefore we can have instability for values of Jo larger 
than 1/4. For R = 2, Ri(y) obtains its minimum value for 
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FIG. 1: The local Richardson number Ri(y) normalized by 
Jo for different values of the parameter R 
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FIG. 2: The potential V(y, c) in equation for c=0 

y = ±oo with Ri(±oo) = 1/4 Jo and therefore instability 
can occur only if Jo < 1. Finally for larger values of 
R, Ri(y) decays exponentially to zero as y becomes large 
and an instability might be present at arbitrarily large 
values of Jo- 

In order to make the second observation we need to 
transform equation (|TJ in to an alternative form. Denot- 
ing E = -k 2 and V(y, c) = U"/(U - c) - J/(U - c) 2 we 
can rewrite the Taylor-Goldstein equation as: 

4>" = [V(y,c) -E]4>. (7) 

Now, if we invert the problem and ask the question, for a 
given real phase velocity c what is the value of fc 2 = —E 
that satisfies the above equation and boundary condi- 
tions, we end up with a Schrodinger problem for a particle 
in a potential well given by V(y,c). Figure El shows the 
resulting singular potential for c = 0, R = 1 and Jo = 2. 
The solution of this particular problem for R = 1 gives 
the stability boundary for the unstable modes given by 
the relation Jo = k(l — k) as shown in [s^. Figure 01 
shows the resulting potential for c = 1 and three values 
of R = 1,2,4. If R < 2 no solution exists that satisfies 
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FIG. 3: The potential V(y, c) in (0 for c=l for three different 
values of R. Only when R > 2 bounded states are allowed. 



the boundary conditions since this problem corresponds 
to finding bounded states in an unbounded potential well 
(V(y, 1) — * — oo for y — > +oo). However if R > 2 we have 
to solve for the eigenstates/energy levels of a particle in 
a finite potential well. For this case a discrete number 
of bounded energy states with "energy levels" E n ex- 
ist. Each state corresponds to a mode with c = 1 and 
wavenumber k n = \J — E n . Since the phase speed of free 
gravity waves is a decreasing function of wavenumber, 
we expect that wavenumbers smaller than y/ —E n will 
have |c| < 1 and could possibly be unstable. Of course, 
whether these solutions indeed represent marginally un- 
stable modes still needs to be demonstrated. Note that 
marginally unstable modes with phase velocity equal to 
the maximum velocity of the shear flow have been re- 
ported before in the literature in interfacial gravity wave 
generation problems [3> HE EH • 

Besides the modes with c = 1, other neutral modes 
that may be marginally unstable should have c in the 
range of U and therefore must be singular. These modes 
can be written in terms of the asymptotic expansion in 
equation ©. For example let <j){y) = A + f a (y) + B + f b (y) 
be the exponentially decaying solution of of the Taylor- 
Goldstein equation ^ f° r V > He (where f a and //, 
are the two Frobenius solutions in [3J|. A + and B + , 
with no loss of generality, are assumed real. The con- 
nection formulas for marginally unstable modes im- 
ply that below the critical layer (y < y c ) the solution 
will be <f)(y) = [A+ cos(s a ir)f a (y) + B + cos(s b n) f b (y)] 
+i[A+ sm(s a ir)f a (y) + B+ sin(s 6 7r)/ 6 (j/)]. However, only 
one linear combination of f a and f b will give an exponen- 
tially decaying solution for y — * — oo. Therefore, if A + 
and B + are such that the real part of cj>(y) is decreasing 
exponentially to zero for y — > — oo the imaginary part of 
4>(y) will increase exponentially. Since both the real and 
the imaginary parts of (j) need to satisfy the boundary 
conditions at infinity it seems unlikely that such a solu- 
tion exists unless one of the two coefficients A + or B + 
is zero. Therefore, a possible marginally unstable mode 
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would be a singular mode that is proportional to only one 
of the two independent Frobenius solutions in equation 
© (i.e. A^ = or _B ± = 0). However the existence of 
such a solution is not guaranteed, so in addition, a suit- 
able value of k 2 and c < 1 needs to be found (if it exists) 
so that both boundary conditions are satisfied. We note 
that the case in which the c = modes form a stability 
boundary corresponds to a special case of this class of 
marginally unstable modes. 



III. NUMERICAL METHODS 

The strategy that we adopt in order to find the unsta- 
ble regions is the following. First we find the location of 
the neutral modes described in the previous section and 
then investigate if these modes are indeed marginally un- 
stable by solving for the growth rate (or for Im{c}) in 
the neighborhood of these modes. Four different numer- 
ical codes were used. Each code integrates the Taylor- 
Goldstein equation using a fourth order Runge-Kutta 
method (tested for different resolutions to verify conver- 
gence) and a shooting method is used to determine the 
eigenvalue in the four different ei genv alue problems that 
we describe in more detail below [37|. 

In the first eigenvalue problem we take c = 0. For a 
given value of Jo we integrate the Taylor- Goldstein equa- 
tion from zero to infinity with initial conditions given by 
one of the two Frobenius solutions 101, and look for the 
eigenvalue k 2 such that the boundary condition at in- 
finity is satisfied. Note that for symmetry reasons, when 
c=0, only the positive y-axis needs to be considered. The 
solution of this problem provides us with the stability 
boundary for the Kelvin-Helmholtz modes, as shown in 
[l7j . We will refer to this problem as eigenvalue prob- 
lem one. The values of Jo that satisfy the condition 
c = for a wavenumber k will be denoted by the curve 
Jo = JKH(k). 

For the second eigenvalue problem we take c = 1 and 
solve the Schrodinger problem for the potential V(y, 1) 
by integrating equation J7J) from — 'oo' to +'oo' and look 
for the eigenvalue E n for which the boundary condition 
at +'oo' is satisfied (where by 'oo' we mean sufficiently 
large value of y). The modes evaluated from this process 
are non-singular neutral modes, and form a curve in the 
(Jo — k) plane that we denote by Jo = Ji(fc). 

In the third problem we are looking for the values of 
c (c / and Im{c} = 0) and k such that one of the 
two independent Frobenius solutions of equation © sat- 
isfies both boundary conditions at infinity. In more de- 
tail the third code begins with an initial "guess" of both 
c and k and integrates from the critical height both for- 
ward and backward up to some sufficiently large posi- 
tive/negative value of y. The initial conditions for cj) and 
cj>' at y = y c ± dy (where dy is a small deviation from 
the critical height) are obtained from equation (© with 
either A ± = or B ± = (both cases were tested) and 
terms up to second order in the expansion are used. As 



argued in the previous section, only the solution that it 
is proportional to one of the two linearly independent so- 
lutions of can satisfy the connection formulas across 
the critical layer (for marginally unstable singular modes) 
and both boundary conditions at infinity. Any other lin- 
ear combination of the Frobenius solutions J3J will lead 
to a singular neutral mode that does not belong to an 
instability boundary. Since we begin with either = 
or = we are only left with two conditions to satisfy 
at ±'oo' and two parameters to change (c and k). With a 
bisection method the code converged to a single pair 
of values of c and k. For one of the two initial conditions 
(say A^ — 0) the code converged to the c = 1 solution, 
and for the second initial condition (say i? ± = 0) the 
code converged to a value of c and k with — f < c < I. 
The solutions with c / 1 comprise a class of singular 
neutral modes forming a curve in the ( Jo — k) plane that 
we will denote by Jo = Jis(k). 

However, the existence of the above mentioned neu- 
tral modes does not guarantee that they form a stability 
boundary. Hence, a code that solves the full eigenvalue 
problem for the complex eigenvalue c was used to ver- 
ify that the above mentioned neutral solutions constitute 
the stability boundaries. This fourth code integrates the 
Taylor Goldstein equation from ±'oo' to zero for some 
given values of Jo and k and a Newton-Raphson method 
is used to find the value of c for which the two solutions 
match at zero. No approximation is used to cross the crit- 
ical height and for this reason this code converged only 
for unstable modes (Im{c} ^ 0) and neutral non-singular 
gravity waves (\c\ > I). Note that this code encounters 
problems with convergence when the examined parame- 
ters are very close to a neutral singular mode stability 
boundary because of the almost singular behavior in the 
neighborhood of the critical height that needs to be re- 
solved. 



IV. RESULTS 
A. Small values of Jo 

We begin with relatively small values of Jo (order unity 
or less) which is the case that has been previously exam- 
ined by 01 an d E3- However in an d 0] the au- 
thors determined the instability boundaries in the ( Jo — k) 
plane by solving the full eigenvalue problem in the exam- 
ined parameter space, finding in that way the regions 
with non-zero growth rates. This is a very difficult task 
since in order to determine the instability region, the en- 
tire (Jo — k) plane needs to be mapped for every value 
of R. Here we use a different approach. We find the 
location of the neutral modes corresponding to the dif- 
ferent eigenvalue problems described in the last section 
and show that these form instability boundaries by solv- 
ing the complex value of the phase velocity only for a few 
values of the global Richardson number. 

In figure 01 we show the results for the case R — 4. As 
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FIG. 4: Panel (a): The stability boundaries for R = 4. In 
region (I) stable gravity waves exist with |c| > 1. In region 
(II) unstable waves with < |Re{c}| < 1 exist. In region (III) 
unstable KH-modes with Re{c} = exist. In region (IV) 
only neutral singular modes exist. The stability boundary 
Ji(fc) is composed of the modes with c — 1. The stability 
boundary Jxff (fc) is composed of the modes with c = 0. The 
stability boundary Jis (k) is composed of the singular neutral 
modes with c / (I. The horizontal lines correspond to the 
values of Jo examined in panels (b,c). Panel (b): The real 
(solid line) and imaginary (dashed line) phase velocity for 
Jo = 0.4. Unstable modes exist for wavenumbers in the range 
ki < k < k\s- Panel (c): The real (solid line) and imaginary 
(dashed line) phase velocity for Jo = 0.2. Unstable modes 
exist for wavenumbers in the range ki < k < kis- In both 
cases Im{c} is approaching zero in the boundaries of the range 
(fci,fcis). All panels share a common :r-axis. 



was found in previous work 

[H 03 the (J - k) plane 
can be divided in four regions (see panel (a)). In the first 
region (I) free gravity waves exist with real phase velocity 
of absolute value greater than one. In the second region 
(II) unstable Holmboe waves exist with the real part of 
the phase velocity smaller than one but different from 
zero. In the third region (III) unstable Kelvin-Hclmholtz 
modes are present with the real part of the phase veloc- 
ity equal to zero. Finally in the fourth region (IV) only 
singular neutral modes exist. The three lines that sepa- 
rate these regions were constructed by finding the neu- 
tral modes with the properties described in the previous 
section. In more detail the curve Ji(k) that separates re- 
gion (I) from region (II) is composed of modes that have 
phase velocity equal to one (c = 1). The curve Jkh (k) 
that separates region (II) from region (III) is composed 
of singular modes with c = 0. Finally the curve Jis(k) 
that separates region (II) from region (IV) is composed 
of singular modes with c ^ 0. 

Further illumination of how the properties of a mode 
change as we increase the wavenumber is gained by look- 
ing at panels (b,c) of figure 0] that show the dispersion 
relation (phase velocity as a function of the wavenum- 
ber) for two values of J (panel (b) J = 0.4 and panel 
(c) Jo = 0.2). A logarithmic scale is used to focus on 
the details of small wave-numbers and small phase ve- 
locities. For the Jo = 0.4 case, for small wavenumbers 
the phase velocity is real and decreases with k until the 
phase velocity becomes equal to one (c = 1 mode) for 
some wavenumber ki ~ 0.05 such that Ji(fci) = 0.4. For 
wavenumbers larger than k± the phase velocity becomes 
complex indicating that the mode with c = 1 belongs 
to the stability boundary. As we further increase the 
wavenumber there is a critical value k\s for which the 
imaginary part of the phase velocity becomes zero again, 
with the real part of the phase velocity remaining finite 
and smaller than one. This mode with wavenumber kis 
can only be a singular neutral mode and can be repre- 
sented by the expansion given in equation with the 
connection formulas given in equation 10}. Indeed the 
difference between the wavenumber where the imaginary 
part of c becomes zero (found by the fourth eigenvalue 
problem) and the wavenumber k\s (found by solving the 
third eigenvalue problem), is found (for this case and all 
the cases examined) to be of the order of the numerical 
error. The small difference is attributed to the difficulty 
in resolving the almost singular region around y c with 
the fourth code when Im{c} is very small. 

For smaller values of the global Richardson number 
(Ri(0) = Jo < 0.25) we have the extra feature that as 
we approach the stability boundary composed of modes 
with c = the real part of the phase velocity decreases 
to zero. It remains zero in region (III) and then starts to 
increase again in region (II) until region (IV) is met. We 
note that in region (IV) a continuum of singular stable 
modes exist and for this reason we do not plot the phase 
velocity after this point. 

The next question we examine is how the stability 
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FIG. 5: The dependence of the instability boundaries on the 
parameter R. The four panels show the stability boundaries 
for four values of R. Panel (a) R = 4, panel (b) R = 2.5, 
panel (c) R = 2.2, panel (d) R — 2. Note that for the last 
case, R = 2, in panel (d) the stability boundary Ji(fe) com- 
posed of modes with the property c = 1 and the boundary 
with singular modes Jis(k) have collapsed together and the 
resulting boundary now defines the region that non-singular 
neutral gravity waves exist. All panels share a common i-axis. 



boundaries change as we vary the value of the param- 
eter R. In figure we plot the stability boundaries for 
four different values of R (panel (a) R — 4, panel (b) 
R = 2.5, panel (c) R = 2.2 and panel (d) R = 2). For 
large values of the parameter R (R > 2) the ( Jo — k) plane 
is divided in four regions determined by the boundaries 
J\, J\s, Jkh as discussed in more detail in figure The 
wavenumbers k of the Holmboe unstable modes are con- 
fined in a finite strip. As we decrease the value of R the 
width of the strip with Holmboe unstable wavenumbers 
becomes smaller (see panel b and c) . When R = 2 (panel 
d) the two boundaries J±s an d Jis collapse together. In 
this case Ji, J\s does not provide a stability boundary but 
continues to separate the region that non-singular gravity 
wave modes exist (on the left of Ji) from the region that 
only singular neutral modes exist with < c < 1 (on the 
right of Ji). The only unstable modes for this case are 
the Kelvin-Helmholtz modes confined in the finite region 
in the left bottom corner. 

We note here that did not find Holmboe instabil- 
ity for values of of R smaller than 2.4. This is probably 
because the width of instability strip becomes very small 
for values of R smaller than 2.4. Furthermore, close to 
the stability boundaries the growth rate is very small 
and as a result it is even harder to find a positive growth 
rate using the method used in 0] when the two stabil- 
ity boundaries are too close. To verify that this strip is 
unstable we plot in figure [5] a close up of the instability 
region and the dispersion relation for Jo = 0.4 for the 
case of R = 2.2. 



B. Higher harmonics 

As we have shown, finding the modes that satisfy c = 1 
can be reduced to finding the eigenstates of a particle in 
a finite potential well. It is well known from quantum 
mechanics that if a potential well is deep enough a finite 
number N > 1 of bound states will exist with N different 
energy levels E n and wavefunctions <j) n (n = 1, 2, . . . , N). 
If the modes with c = 1 form one of the stability bound- 
aries for the Holmboe instability, it is natural to ask what 
happens for the case that Jo is large enough (i.e. the 
potential well in the Schrodinger problem Q is deep 
enough) so that more than one mode satisfies the con- 
dition c = 1 ( i.e. more than one bounded eigenstate is 
allowed). H 

To investigate this we solve for the modes with c = 1 
for values of Jo up to 10 2 . We find that the higher har- 
monics of the gravity waves with c = 1 indeed exist and 
correspond to new stability boundaries that we denote as 
J n (k). The second harmonic (n — 2) with c — 1 appears 
for Jo ~ 28 and the third harmonic (n = 3) for Jo ~ 80. 
We also find that for each c = 1 mode a marginally sin- 
gular neutral mode with < c < 1 also exists. These 
modes form new curves in the ( Jq — k) plane given by 
Jns(k). These curves J n (k) and J n s, define new stripes 
of unstable regions. 
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FIG. 6: The stability boundaries panel (a) and the dispersion 
relation (Jo = 0.4) panel (b) for the case that R — 2.2. In 
panel (b) solid line gives the real part of c and the dashed line 
gives the imaginary part of c. Unstable Holmboe modes exist 
for this case also but are confined in a very narrow region. 




FIG. 7: Instability diagram for larger values of Jo and R = 4. 
Solid lines show J„(k), dashed lines show J„s(k). The stable 
regions are marked with S. Regions marked with U n are 
regions that the n-th harmonic becomes unstable. 
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FIG. 8: Panel (a):The stability boundaries for the second 
harmonic for R = 4. Solid line gives J2(k), dashed line gives 
Jzs(fe)- Panel (b): The real part of the phase velocity (solid 
line) and the imaginary part of the phase velocity (dashed 
line) for Jo = 35. 



modes exist in the regions between J n (k) and J n s{k) by 
solving for the complex phase velocity of the modes in 
these regions. In figures [S] and [5] (panel (a)) we show 
a close up of the stability boundaries J2{k), Jzsik) an d 
Ja(k), J3s(k) respectively, and the dispersion relation in 
panel (b) . Figure 00 and [5] demonstrate how the second 
and the third gravity wave harmonic become unstable. 



In figure [7| we plot the stability domain for values of 
Jo up 100 for R = 4. The stable regions are marked with 
S and the regions that the n-th harmonic becomes un- 
stable are marked with U n . The solid lines show J n (k) 
and correspond to modes with c = 1 and the dashed 
lines show J n s{k) and correspond to marginally unsta- 
ble singular modes. As before we verify that unstable 



We note finally that for a fixed value of Jo the high- 
est harmonic has the largest growth rate. In figure I1UI 
we show the growth rate 7 = fclmjc} as a function of 
the wavenumber for R = 4 and two different values of 
Jo. The growth rate is significantly larger for the high- 
est harmonic. This makes the higher unstable harmon- 
ics possible to be detected experimentally, since it is the 
most unstable modes that are usually observed. 
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FIG. 9: Panel (a). The stability boundaries for the third 
harmonic for R — 4. Solid line shows Js(fc) and dashed line 
shows J3s(k). Panel(b) The real part of the phase veloc- 
ity (solid line) and the Imaginary part of the phase velocity 
(dashed line) for Jo = 85. 



V. CONCLUSIONS AND PHYSICAL 
DESCRIPTION 

In this paper we investigated Holmboe's instability in 
stratified shear flows for smooth density and velocity pro- 
files. Using a specific model of the density and velocity 
profiles |l7lll8j |. we were able to determine the instability 
regions for large values of the global Richardson number 
including regions not previously noted in the literature. 

Focusing on the nature of the stability boundaries that 
enclose the Holmboe unstable modes, we have shown that 
for moderate values of Jo and when the density stratifi- 
cation length scale is sufficiently smaller than the shear 
length scale (i.e. R > 2) the ( J — k) plane is divided into 
four regions: (I) a region where neutral gravity waves ex- 
ist (i.e. modes with real phase velocity outside the range 
of the velocity shear); (II) a region where unstable trav- 
eling waves exist (i.e. Holmboe modes); (III) a region 
where unstable modes with the real part of the phase ve- 




Tb=3 



= 45 J 



rv= 1 



(b) J =15 



FIG. 10: The growth rate for the different unstable harmon- 
ics. Panel (a) Jo = 45, R = 4. Panel(b) Jo = 15, R = 4. For 
the smallest harmonics for which the Im{c} was very small the 
code was able to find the growth rate only around the peak 
and not close to the stability boundaries. For the Jo = 45 case 
in panel (a) the growth rate of the first harmonic (n = 1) was 
too small for the code to be able to resolve the critical height 
and the dark line just indicates the location of the instability 
region based on the stability boundaries shown in figured 



locity being zero exist (i.e. Kelvin Helmholtz modes); 
(IV) and finally a region where only singular neutral 
modes exist. We determined the modes that comprise 
the boundaries and separate these four classes of modes. 
For small Jo (Jo < 1/4) Kelvin-Helmholtz modes exist. 
They are restricted in a bounded domain enclosed by the 
boundary Jo = JKn(k) determined by the modes with 
zero phase velocity. For larger Jo , a strip of Holmboe un- 
stable modes exist. This instability strip is determined 
by the two boundaries Jo = Ji(fc) and Jo = Jis(fc). In 
more detail, we have demonstrated that for a given value 
of the global Richardson number Jo the unstable modes 
have wavenumbers k that lie in the range k\ < k < kis 
where k\ is the wavenumber of the mode that has phase 
velocity equal to one and satisfies Jo = Ji(fci) and his is 
the wavenumber of a singular marginally unstable neu- 
tral mode. k±s satisfies Jo = Jis(kis)- We have shown 
how the value of ki and his can be determined to desired 
accuracy (and therefore how to construct the boundaries 
Ji(fc) and Jis(fc)) by solving a Schrodinger eigenvalue 
problem for a particle in a potential well. 

For large values of the global Richardson number more 
than one strip of instability may appear. Each new 
instability strip is confined between the pair of curves 
Jo = J n (k) and J = J n s(k). As before J n {k) can be de- 
termined for a given value of Jo by finding the wavenum- 
ber k n of the n— th harmonic of the gravity wave spec- 
trum with c = 1 and J n s(k) is determined by finding 
the wavenumber k n s of the singular marginally unstable 
mode that is met first as we increase the wavenumber 
from k n . 

The results of this paper show a deep connection be- 
tween the free gravity wave spectrum and the Holmboe 
unstable waves. The overall picture looks as follows. For 
large stratification and small wavenumbers the free grav- 
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ity wave spectrum is only slightly modified by the shear 
and is composed of a discrete finite number N (N > 1) of 
stable modes (0 n , c„) with n = 1, 2, . . . , N. As the wave 
number is increased the phase speed decreases roughly as 
the square root of the wave number. As the phase speed 
approaches the maximum shear velocity (sup{t/}), the 
functional form of the stratification plays a crucial role. 
If the stratification length scale is larger than a critical 
value (more precisely if lim^oo Ri(y) — > oo) the phase 
speed approaches sup{J7} asymptotically as k — > oo, al- 
ways remaining outside the range of U. In this case 
the only unstable wavenumbers are the Kelvin- Helmholtz 
modes that have zero phase velocity and appear only for a 
finite range of the global Richardson number. If the strat- 
ification length scale is smaller than this critical value 
(more precisely if lim^oo Ri(y) — ► 0), the phase speed 
of the waves reaches the value of the maximum wind 
speed for a finite value of the wavenumber k = k n . As 
we increase the wavenumber further, the modes become 
unstable with the real part of the phase velocity smaller 
than one. The instability persists up to another critical 
value of the wavenumber k — k n s ■ The mode with this 
wavenumber exhibits a singular behavior at the critical 
height. For wavenumbers smaller than k n s a continuum 
of singular neutral modes exist. 

The findings in this work suggest several possibilities 
for further investigation. An important aspect that needs 
to be investigated is the effect of viscosity and diffusiv- 
ity. As the work in [T^. |20| has shown, the presence of 
finite viscosity and diffusivity decreases the growth rate 
and this effect is expected to be stronger for the higher 



harmonics. This may limit attempts to detect the higher 
harmonic modes via numerical and experimental stud- 
ies. Therefore, an investigation of the viscosity effects is 
further needed to find the critical Reynolds number for 
which the higher modes become unstable and this should 
be the next step in the study of the higher harmonics in- 
stability. Nonetheless, we note that in most astrophysical 
and geophysical flows, where Holmboe's instability is of 
great importance, the Reynolds and Peclet numbers are 
large enough that the inviscid flow examined here is a 
good approximation, and depending on the stratification 
the higher harmonic modes might dominate. 

Finally, besides viscosity and dissipation, three dimen- 
sional effects, nonlinear evolution and saturation, as well 
as physical interpretation of these results through sim- 
plified layer models are all important issues that we plan 
to investigate in our future work. 
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